; docformat = 'rst'
;
;+
;
; :Purpose:
;   Calculate mid-point hybrid-coordinate pressure fields.
;   For the ith level:
;   
;     P_i=A_i*P_0 + B_i*P_surface
;   
;   If /half is set, it is assumed that A (and B) define the coordinate system at
;   grid cell boundaries. Mid-point pressure values will the be output, using:
;     
;     P_i=(P_{i-1/2} + P_{i+1/2})/2
;
; :Inputs:
;   A: (input, floating, required), hybrid A values. If P0 is set, A is dimensionless.
;     If A P0 is not set, A must be in Pa.
;   B: (input, floating, optional), hybrid B values. Dimensionless. 
;     If left out, surface pressure independent (fixed) pressure levels are output. 
;
; :Keywords:
;   P0: (input, optional, floating) Reference pressure to multipy the A coordinate
;   Half: (input, optional, boolean) Set to True if A and B define coordinates at grid cell
;     boundaries. The output will still be at grid cell centres.
;     Note that if true, the output will have one fewer vertical elements than A and B.
;
; :Outputs:
;   3D pressure field (lon, lat, level)
;
; :Example::
; 
;   IDL> Pressure=mr_hybrid_coords(A, B, PS, P0=1.e5)
;
; :History:
; 	Written by: Matt Rigby, University of Bristol, Jul 16, 2013
;
;-
Function mr_hybrid_coords, A, B, PS, P0=P0, half=half

  on_error, 2


  ;Define dimensions and output array  
  LevSize=n_elements(A)
  LonSize=n_elements(PS[*, 0])
  LatSize=n_elements(PS[0, *])
  P=dblarr(LonSize, LatSize, LevSize)

  ;Check default inputs
  if n_elements(B) gt 0 then begin
    if n_elements(B) ne LevSize then message, "A and B must have same dimensions"
    has_B=1
  endif else begin
    has_B=0
  endelse
  
  ;Check if A is dimensionless
  if ~keyword_set(P0) then begin
    P0=1.
  endif
  
  ;Calculate pressure
  For LevI=0, LevSize-1 do begin
    P[*, *, LevI]=A[LevI]*P0
    if has_B then P[*, *, LevI] += B[LevI]*PS
  endfor

  ;If /half, calculate mid-point pressure levels
  if keyword_set(half) then begin
    P_out=dblarr(LonSize, LatSize, LevSize-1)
    for LevI=0, LevSize-2 do P_out[*, *, LevI]=0.5d*(P[*, *, LevI] + P[*, *, LevI+1])
    P=temporary(P_out)
  endif

  return, P

end